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We present Quantum Monte Carlo calculations with up to A'^ = 576 000 interacting bosons in 
a quasi two-dimensional trap geometry closely related to recent experiments with atomic gases. 
The density profile of the gas and the non-classical moment of inertia yield intrinsic signatures for 
the Kosterlitz-Thouless transition temperature Tkt- From the reduced one-body density matrix, 
we compute the condensate fraction, which is quite large for small systems. It decreases slowly 
with increasing system sizes, vanishing in the thermodynamic limit. We interpret our data in the 
framework of the local-density approximation, and point out the relevance of our results for the 
analysis of experiments. 
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Phase transitions in two-dimensional systems with a 
continuous order parameter are of special interest be- 
cause long-range order is absent at finite temperatures, 
as a consequence of the Mermin- Wagner theorem fl|, Q ■ 
Instead, a Kosterlitz-Thouless phase transition [3| can 
separate the high-temperature disordered phase with ex- 
ponential decay of the order parameter from a low- 
temperature ordered phase with algebraically decaying 
order parameter, which was proposed by Berezinskii 

Recently, the Kosterlitz-Thouless phase transition was 
observed in trapped quasi two-dimensional Bose gases of 
^^Rb atoms 0, 0] , but the interpretation of the experi- 
mental data was rendered difficult because of relatively 
strong interactions in two dimensions, the trap confine- 
ment, and pronounced finite-size corrections. Moreover, 
the ideal two-dimensional trapped Bose gas is not a good 
vantage point to approach the weakly interacting gas: the 
former Bose-condenses at finite temperature (with a di- 
verging density in the center), whereas the latter becomes 
a superfluid, but with vanishing condensate fraction. 

Quantum Monte Carlo (QMC) calculations allow us 
to compute the thermodynamic properties of interacting 
Bose systems Q for a finite number N of particles, and 
for a wide range of microscopic interaction parameters. 
The calculations are practically free of systematic errors. 
For trapped atomic gases, very large particle numbers can 
be handled, in confining geometries relevant to current 
experiments 0, [^, . The ideal Bose gas enters these 
calculations in an exact way, and only pair interactions 
give rise to the usual Metropolis rejection process 

In this paper, we present QMC calculations of three- 
dimensional trapped bosons interacting with an s-wavc 
pseudopotential in a pancake-shaped harmonic trap with 
frequencies uj = uj^. ^ <C lOz- Our particle numbers 
range from N = 2000 to TV = 576 000, the latter exceed- 



ing current experiments by more than one order of magni- 
tude. The diagonal many-particle density matrix directly 
yields the density profile and the non-classical moment of 
inertia. Both allow us to locate the phase transition. We 
also determine the condensate fraction explicitly from the 
largest eigenvalue of the reduced off-diagonal one-particle 
density matrix. In a trap, this calculation is more com- 
plicated than in a homogeneous system, where the fc = 
groundstate wavefunction is trivially known and where 
the groundstate occupation governs the long-range be- 
havior of the off-diagonal one-particle density matrix. 

The trapped ideal two-dimensional Bose gas shows a 
Bose-Einstein transition [l^ at a non-vanishing tem- 
perature Tbec ^ ■\/ QNlo'^ /tt"^ (we choose units with 
ft = m = w = 1, where m is the atomic mass). The cen- 
tral density diverges logarithmically with N (see Fig.[T]). 
This implies that interaction effects play a much more 
pronounced role, even above the transition temperature, 
than in three-dimensional traps. 

As in the experiment , we allow for a finite extension 
LOz of the trap in the z-direction, keeping the level spac- 
ing on the order of the temperature: lOz = 0.55 Tg^^. 
In the many-body density matrix, the z-dependence is 
dominated by a (normalized) single-particle contribu- 
tion, p(z, z'), which separates out, and the effective two- 
dimensional interaction strength is given by 



9 = 



dz [p(z,z)]^ 



(1) 



where gq is the three-dimensional s-wave scattering 
length. For particles distributed in z according to the 
groundstate of the harmonic oscillator, Eq. ([T]) reduces 
to g = g = ao-v/SmjJ. For our simulations, we have used 
the experimental value 5 = 0.13 Q, however, the actual 
value of g can be obtained directly from the computed 
density profile in the z-direction. 



2 



r=o 

QMC 
Boltzin. 




n,.A 



T=0 
QMC 
Boltzm. 




FIG. 1: Two-dimensional density profile n{r)X^ at T = Tbec 
for A'^ = 576 000, compared to the saturation density of the 
ideal Bose gas and the density profile of an ideal gas of dis- 
tinguishable particles. The inset compares the density profile 
in z to the groundstate distribution of the harmonic oscillator 
and to the ideal gas of distinguishable particles. 



We study the anisotropic trap at temperatures compa- 
rable to TggQ where it is indeed quasi two-dimensional 
because the extension in z is comparable to the de- 
Broghe wavelength A = yJ^-KjraT. Formally, the 
three-dimensional Bose-Einstein transition temperature 

ever, the three-dimensional limit requires that uizjio re- 
mains constant independent of the system size in contrast 
to our quasi two-dimensional limit where uJz/^ ^ N^^^. 

Recent numerical calculations have determined the 
critical density ric at the Kosterlitz-Thouless transition 
in the weakly interacting two-dimensional homogeneous 
Bose gas of density n, 
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The interaction g enters this expression only logarithmi- 
cally, and the differences between the actual g and g, of 
the order of 40% at Tkt, only results in a 6% shift in the 
critical density. 

In the trapped Bose gas, within the local-density ap- 
proximation, the transition takes place when the central 
density n(0) equals the critical density of the homoge- 
neous gas, in our case nc(0)A^ ~ 8. Mean-field theory 
14| predicts that the Kosterlitz-Thouless transition is 



somewhat below TggQ. 
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The mean- field value of Tkt, together with the numerical 
value of Eq. ([2]) for the critical density in the center of 
the trap allows to determine the critical temperature of 
the Kosterlitz-Thouless transition in the trap, in our case 



FIG. 2: Two-dimensional density profile n(r)A^ at tempera- 
ture T = O.ST^EC for N = 576 000 (thick line), compared 
to the Thomas-Fermi profile of Eq. Q (with g = 0.109, 
dashed line). The ansatz of Eq. ([5]) for the superfluid den- 
sity Ps{r), with the universal jump at r = Vc, corresponds to 
the shaded region. The inset compares the density profile in 
the tightly confined z-direction to the groundstate distribu- 
tion of the harmonic oscillator and the distribution of an ideal 
gas of distinguishable particles. 



In Fig. [U we show the two-dimensional density pro- 
file n{r) with r ~ \J -V'ip' from our QMC calculations 
at T = T^EC for iV = 576 000. We also illustrate the 
large deviations from the saturation density of the two- 
dimensional ideal Bose gas [ll| . Indeed, the density pro- 
file is closer to that of ideal quantum Boltzmann particles 
described by the density matrix of the harmonic oscilla- 
tor. The density is everywhere below the critical value, 
confirming that the interacting gas remains in its high- 
temperature phase at lower temperatures than the ideal 
Bose gas. 

In Fig. [21 we show the analogous density profile at tem- 
perature T/TegQ = 0.5, again for N = 576 000 particles. 
The central density is now well in excess of the criti- 
cal value of Eq. ([2]). We may define a "critical radius" 
Tc, which separates the "inner region" of the trap, with 
r < Tc and n{r) > ric, from an "outer region" with r > 7'c 
and with n{r) < ric- In the local-density approximation, 
the inner region is in the superfiuid phase, whereas the 
outer region is normal. At the critical radius, the density 
is at the Kosterlitz-Thouless phase-transition tempera- 
ture. In the inner region, n{r) is very well described by 
a Thomas-Fermi profile 



n[r) = nifi) — 



2.9 



(4) 



0.75 TOgc- 



with the effective two-dimensional interaction parameter 
at this temperature g = 0.107, obtained, via Eq. ([1]), di- 
rectly from the density profile in z (see the inset of Fig.[2]) . 
The latter is wider than the groundstate distribution of 
the harmonic oscillator, so that g is smaller than g. The 
density profile in r, whose width depends linearly on g~^ , 
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FIG. 3: Left: Densities ricX^ and n(0)A^ vs. T/T^ec for 
the quasi two-dimensional gas with = 576 000. The in- 
tersection ol both curves leads to a transition temperature 
Tkt — 0.70. Right: central curvature k compared to the 
groundstate Thomas-Fermi curvature, Eq. (U), with g = g 
and with g corresponding to an ideal gas of distinguishable 
particles. The central curvature changes slope at Tkt. 



is more sensitive to the detailed value of the interaction 
than the transition temperature, which decreases with 
the logarithm of g. 

In Fig. 131 we plot the central density n(0)A^ and also 
the (central) curvature k = —{X'^/P)dn{r)/d{r^)\r=Q- 
The curvature of the Thomas-Fermi profile (in Eq. (|4|)) 
is K = uP'-Kjg with g = g a.t very low temperature. The 
curvature increases (the profile becomes narrower) with 
T because particles spread out farther in the z-direction. 
Above the critical temperature, however, the curvature 
decreases (the profile becomes wider), as is natural for a 
thermal gas, with k cx n(0)A^. The curvature plot pro- 
vides an intrinsic signature of the phase transition, at a 
temperature T/T^^^ ~ 0.70, which agrees nicely with 
the temperature at which the central density passes the 
critical value Eq. ([2]). We have also studied smaller sys- 
tems (with N = 2250,9000,36 000, and 144 000) at un- 
changed values of T/T^^^ and Wz/Tg^^, but found only 
very small variations in the density profiles. Our value 
for the critical temperature, Tkt — 0.70 TgEQ is close to 
the mean- field formula of Eq. ([3]) , and somewhat higher 
than the experimental value 0] ~ 0.46(3) TggQ. 

The low-temperature phase below the Kosterlitz- 
Thouless transition is a superfluid. For a homogeneous 
system, the superfluid fraction can be probed through the 
response to boundary conditions, and easily computed 
within the path-integral formalism, through the winding- 
number formula Q. Likewise, a trapped superfluid does 
not respond to an infinitely slow rotation of a trap lead- 
ing to a non-classical moment of inertia, /nc, which is 
smaller than the classical value Ic\ = J drr^n(r). The 
non-classical moment of inertia can again be computed 
from the diagonal elements of the density matrix [15| . In 
a homogeneous system, the ratio of the non-classical mo- 



0.7 



-d 

o 



I 



0.4 0.6 



0.8 



BEC 



0.4 0.6 



0.8 



FIG. 4: Left: Non-classical moment of inertia /nc//ci vs. 
T/Tisc for N = 9000 (crosses) and = 144 000 (stars) com- 
pared to the ansatz of Eq. ([5]) (squares) . Right: Condensate 
fraction for particle numbers ranging from A'^ — 2250 (crosses) 
to iV = 144 000 (squares). 



mcnt to the classical moment equals the normal fraction. 
In Fig. m we show that a superfluid phase emerges below 
T ~ 0.70 TggQ, and that Inc/Ici remains different from 
unity, independent on system size. 

To interpret our data for the non-classical moment of 
inertia, we observe that in an infinite homogeneous sys- 
tem, at the Kosterlitz-Thouless transition, the superfluid 
density develops a universal jump[l6|. Aps = 2toTkt/7'', 
and the superfluid mass and the moment of inertia are 
both discontinuous. In the trap, the spatial structure 
smears out these discontinuities |14l | , but in local-density 
approximation, as mentioned, the gas is critical at the 
critical radius Vc. Therefore, we expect a normal phase 
beyond rc, and a superfluid for r < Tc, with a jump 
of the superfluid density taking place at this radius and 
the superfluid density vanishing for r > Tc. For our pa- 
rameters, the superfluid fraction at the critical radius 
is Ps{rc)/n{rc) = 2mT/n(rc)'!T ~ 0.5. We can continue 
the superfluid density Ps{r) into the inner region by a 
Thomas-Fermi profile: 



(1 - r^/ri) /2g + 2mT/7r, for r < 



Psir) 



2niT/TT for r 
for r > Vc 



(5) 

(see Fig. The non-classical moment of inertia I^c = 
/ dr r^[n{r) ~ Ps('')], computed using Eq. (O and the 
computed density profile n(r), agrees excellently with our 
data (see Fig. [4|). 

In the low-temperature phase of a two-dimensional su- 
perfluid, the condensate fraction vanishes in the thermo- 
dynamic limit. In a homogeneous system, the ground- 
state has zero momentum, and the fraction of parti- 
cles occupying this state can be computed from the 
long-distance behavior of the non-diagonal one-body 
density matrix, p^^^ (r, r'; /3) (see 0). In an inhomo- 
geneous system, the groundstate eigenfunction of the 
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one-body density matrix is no longer completely deter- 
mined by symmetry. Still, in the rotationally symmet- 
ric trap the one-body density matrix is block-diagonal 
with respect to the Fourier components / of the angle 
between r and r'. Projection onto the Fourier com- 
ponents yields onc-dimensional matrices, r'; 
which can be discretized more easily than the bigger ma- 
trix p(^)(r, r'; /?). The condensate fraction, Nq/N, cor- 
responds to the largest eigenvalue with I ~ 0. Con- 
densate wavefunctions computed this way for the three- 
dimensional trapped Bose gas, with Legendre polynomi- 
als replacing the Fourier components, closely agree with 
the solution of the three-dimensional Gross-Pitaeveskii 
equation (see Figure S] shows that the conden- 

sate fraction of our quasi two-dimensional system is 
rather large, but it decays algebraically with system 
size: Nq/N - j^-v{t)/2^ r^^ie exponent -qiT) depends 
on the temperature; we obtain 77(0.70 Tgg(^) k, 0.5 and 
77(0.67 rgg^-,) w 0.2. Precisely at the critical temperature, 
we expect ri{TKT) — 1/4 which implies that the critical 
temperature is between 0.67 Tg^^ and 0.70 Tg^j-,, com- 
parable with our previous estimate, Tkt — 0.70 Tg^^, 
based on the occurence of a non-classical moment of in- 
ertia. 

Our results are in qualitative agreement with recent ex- 
periments 0. However, they found a lower critical tem- 
perature = 0.46(3) TggQ, and an almost Gaussian 
shape of the density profile above Tkt- In the specific 
experiment, the presence of several planes in the optical 
lattice renders the estimation of the critical temperature 
and the particle numbers difficult. More generally, we 
see from Fig.[Tl even in the normal state away from Tkt, 
that the density profile deviates from the classical Gaus- 
sian distribution as soon as 7^(^)A^ > 1. Therefore, only 
a small part of the density distribution can be used to 
gauge the temperature close to Tkt- Fitting a larger 
part of the distribution to a Gaussian typically yields a 
smaller width than the classical thermal distribution and 
therefore underestimates the true temperature of the sys- 
tem. This problem is even more pronounced when one 
analyzes column densities rather than radial profiles. 

Notably, in the trap, the universal jump in the super- 
fluid density of a two-dimensional superfluid at the criti- 
cal temperature does not induce a significant discontinu- 
ity in the inertial response, as in two-dimensional films 
[13]. Nevertheless, the universal jump determines Ps{r) 
at the radius r = Tc where the gas is locally critical. 
Based on the local-density approximation, we proposed 
a superfluid density proflle, Eq. ([5]), which continues the 
Psir) from r = Tc into the superfluid inner region. It 
depends on a single parameter Tc whose value can be 
determined directly from the density profile. The non- 
classical moment of inertia calculated from the superfiuid 
density profile, Eq. ([5|), is in excellent agreement with a 
direct computation of this quantity. Further, it is re- 



markable that in the finite Kosterlitz-Thouless system, 
the condensate fraction, which must vanish for an infinite 
system, is still rather large, even close to the transition 
temperature. The fact that the groundstate wavefunc- 
tion of size ^ Tc remains macroscopically occupied for 
systems with particle number < 10^ implies that the 
coherence of the atoms is neither destroyed by interparti- 
cle interactions nor by fluctuations, essential for building 
continuous and coherent sources of matter waves in lower 
dimensions [l9| . 
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